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We present a numerical scheme for tiie Hubbard model that throws light on the rather esoteric 
nature of the Upper and Lower Hubbard bands that have been invoked often in literature. We 
present a self consistent solution of the ladder diagram equations for the Hubbard model, and show 
that these provide, at least in the limit of low densities of particles, a vivid picture of the Hubbard 
split bands. We also address the currently topical problem of decay of the doublon states that 
' • are measured in optical trap studies, using the ladder scheme and also by an exact two particle 

' ' calculation of a relevant Greens function. 

o : 

^ ■ I. MOTIVATION AND INTRODUCTION 

X) : 

' Hubbard's introduction of split bands in Ref.Q, i.e. the so called upper Hubbard band (UHB) and the lower 
Hubbard band (LHB), is one of the most important qualitative ideas in the theory of correlated electrons. Their 
t ' origin is the idea that since the energy levels of the atomic limit show two sets of states, one at a; ~ and another at 
w ~ J7 as in Eq. below, the formation of a crystal would broaden these levels into two sets of sub-bands. These 
' sub-bands were originally discussed by Hubbard using a non perturbative technique, that has the advantage of being 
I , exact in the limit of vanishing bandwidth W 0, i.e. the atomic limit. However, the technique failed to produce 
^ ■ a Fermi liquid for weak couplings, as one expects physically. This failure led to severe early criticism of Hubbard's 
^ , work?-. The problem of reconciling Fermi liquids with the local picture developed by Hubbard, leading to the split 
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bands, is of great importance in the physics of strong correlations. The one exception is the dynamical mean field 
theory that gives a good account of the sub-band formation, especially in the proximity of half filling^>^. However, 
G , away from half filling, the picture is obscure and remains largely unresolved. It is this task that we address in the 
■ present work. We study the ladder diagrams that are argued to be exact at low densities, sharpen the argument 
, for their validity in terms of the self energy, and show that at least in this limit, the concept of the split bands is 
O ' completely consistent with the Fermi liquid picture. The numerical solution of the ladder diagrams is carried out in 
' a self consistent way and shows the emergence of the Hubbard split bands for large enough U/W. These merge for 
weak couplings and our results give a vivid picture of the crossover from weak to intermediate to strong coupling. 
I The self energy is momentum and also frequency dependent in the ladder scheme, and for low densities provides a 

^ full picture of the renormalization processes that occur at arbitrarily large interaction scale U. In particular we see 
that the spectral function shows a low lying feature and a high energy ^ 0{U) feature, with spectral weights that 
are equal to 1 — ^ and ^ respectively. It is seen that every single added particle thus depletes the weight of the LHB 
' and adds to the UHB, thereby accomplishing a "long range spectral transfer"- that has been described in literature 
as "Mottness"-^!^. 

. The momentum space occupancy m{k) = (c^^Cfco-) is computed and it is usefully broken up into three parts Eq. (|16l) . 
The occupied part mi(fc) in Eq. corresponding to occupied states that are automatically inside the LHB, the 

unoccupied LHB part m2(fc) corresponding to unoccupied LHB states, and the unoccupied UHB part m3(fc). In the 
limit oiU oo, only uii and TO2 survive, and this projection gives an exact view of the physics of the t-J model as well 
in the low density limit. At low densities we find that the ladder diagrams lead to a Luttinger Ward compliant Fermi 
surface, and this Fermi surface survives the limit [/ — > oo. Thus even in this limit of extreme correlations J7 — > oo, 
adiabatic continuity to the Fermi gas holds. Therefore we have a useful and concrete alternative to the extreme 
Ci ' coupling ideas proposed in work by one us^, where a different Fermi volume emerges at all densities, including the 
lowest ones. 

One contemporary context for the Hubbard split bands is the problem of high Tc superconductors, here Anderson^ 
has eloquently argued that for large U, one can confine attention to carriers in the LHB, with the UHB pushed out 
of the range of relevant states. Given this projection to the LHB, the charge carriers inherit exotic properties such as 
spin charge separation, and also a new interaction, namely the super exchange that comes with a scale of jU . We 
see that at least at low densities where the ladder scheme is valid, the LHB does separate out cleanly for U > W, but 
the carriers are yet subject to Fermi liquid behaviour. 

Another recent context for motivating this work is the study of the Hubbard model far away from equilibrium with 
cold atom realization^!^, where the carriers in the UHB are optically excited, and their lifetime studied by measuring 
the overlap of the excited state with the initial state. We find that a calculation of a related correlation function is 
possible in the Fermi liquid at low densities, albeit in a close to equilibrium situation unlike the experiments. We are 
also able to exhibit the correlation function exactly for a pair of particles in the Hubbard band. Interestingly, the 
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resulting life times show some similarity in functional dependence to those found in experiment, although with a very 
different time scale. 



II. LADDER SCHEME EQUATIONS AT LOW DENSITY 

The ladder scheme for the Greens function Ref s . [l0l- [T2lfl3fTH corresponds to convoluting a particle-particle ladder 
scattering amplitude T{Q) with a single Greens function G{k) to form the self energy as follows: 

r(Q) = ^ 



i + c/n(Q)' 



^(^) = ^EG(p)r(fc + p)- (1) 



Here Ng, are the number of sites and electrons, n — N^/Ng is the electron number density, and we use the notation 
k = {k,iujk) with imaginary odd frequencies uJk = n^{2k + 1) of the finite temperature field theoryi^ for Fermions, 

and reserve the capital letters for Bosonic frequencies, e.g. Q — {Q, iQ^) and fl^, = Here the summation over p 

represents a sum over the vector component and also the imaginary frequency. A paramagnetic state is assumed and 
the spin label is suppressed for brevity. In addition to Eq. ([T]), we have the Dyson equation G~^{k) — G'^^(fc) — 
with the usual non interacting Greens function G(^^(fc) = iuik — + /i. Thus the ladder scheme is a self consistent 
non linear scheme that needs to be solved numerically for the various objects G{k), T,{k), Il{K). We can solve for 
the Dyson equation in the ladder scheme iteratively: 



For example in the first step we can calculate the scattering amplitude (and self energy) using Go and use Dyson's 
eqn to obtain a new Green's function we call Gi: 

G^\k) = G,\k)~-l-^Goip)-—^r^^-^f-—— -. (3) 

P^'^p ^ + jw:J2qGo{q)Go{k+p-q) 

We may continue and compute G2(fc) using Gi(fc) to recompute the self energy (i.e. the second term in Eq. ([3])), 
and repeat this process iteratively to obtain G(fc) = lim„_j>oo G„(fc). The difference between Gi(fc) and the fully self 
consistent G(fc) arises from the repeated renormalizations implicit in the full equations, and this brings about the self 
consistent broadening of several sharp features that arise in Gi(fc). In Fig. ([3]) we discuss the difference in the spectral 
functions from these two theories as an illustration of this phenomenon. 

Alternatively we start by introducing spectral representations for the various quantities of physical interes t^^i"*^^ : 



G(fc,jcjfc) ^ j dv 



PG{k,iy) 
ij(k,iuJk) = U — -\- I dv 



2 J iuik — V 

r(Q,zf7Q) = C/+ / dvP^^. (4) 
J i\Iq-v 

The spectral functions Pt{Q, v) etc have a compact support and are therefore convenient for numerical integration on 
a suitably discretized grid of frequencies. The numerical solution is performed after using a spectral representation 
for various physical quantities. We first turn the Dyson equation Eq. ([2]) into a non linear integral equation for the 
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spectral function from Eq. (|4]) as follows: 

* p 

MQ.n) ^ -Sfsm — ,5, 

(1 + ufle n(Q, n))" + (»t;pn(i?,si))2 

with f{uj) and nsito) as the Fermi and Bosc distribution functions [exp ± and Re Il{Q,il.) defined as the 

Hilbert transform of pYi{Q,v), i-e. 



Re n(Q, f7) = -P y diy 



A. Low density limit and self energy sum rule 

The original argument for the ladder schemeiSili is that it is exact in the low density limit. This argument is 
borrowed from the theory of nuclear matter, where Brueckner— originally argued that at any order n of perturbation 
theory for the ground state energy (i.e. the Goldstone diagrams), the dominant diagrams are those with the smallest 
number of downward lines of holes. Topologically there need to be at least two such hole lines in the free energy 
diagrams. The particle particle ladder diagrams have only two hole lines at any order. Thus the ladder diagrams 
dominate all others at each order in perturbation theory. Importantly for nuclear matter, this logic shows that the 
large (divergent) two body interaction is not a problem, it is cut off by these ladders, giving in the end an expansion 
in a dimensionless parameter obtained by combining the two body scattering length with the average inter particle 
separation. A parallel argument for bosons was provided by Lee, Huang and YangiS.. The Kanamori- Galitskii papers 
implement this idea for the Feynman diagrams, where one has additionally hole hole scattering, in addition to particle 
particle ladders- for structural reasons that distinguish the Ferynman diagrams from the Goldstone ones. However 
these extra terms do not detract from the particle particle ladders that cohabit the Feynman series and provide a 
particular 0{n^) correction term. 

The reader would note that the above argument is rather indirect, in particular it gives us no clue to why we should 
accept the self energy that emerges from this scheme as exact. In this context, it is useful to note that the self energy 
satisfies an exact series of sum rules^ i^^'^^ , of which the lowest is 

so(fc)= J dup^{k,u)^U^ X (6) 

where the RHS is independent of k. Note that this sum rule is valid for arbitrarily large U and at all densities. We 
can use this as a check of our calculation by testing for the k independence of the computed LHS, and also monitor 
its weight relative to the RHS. The self consistent solution of the ladder diagrams contain the low density limit and 
also provide some uncontrolled results at higher densities, and it is important to know the limit on density to which 
we can trust these results. Fig. ([TJ gives details of this test for the ladder diagrams. For higher densities the ladder 
diagram theory is systematically wrong for the O(n^) term, since we can show analytically that at large U and low 
density so(k) = C/^ x + 0(U) in contrast to Eq. 



B. The Atomic Limit 



We discuss briefly the atomic limit, i.e. a limit where U remains finite but the band width W ^ Q, this is the limit 
where one can solve for the Greens function exactly quite simply. 

^Atormc = U X - + X ^ ^ (7) 

2 tuj + ^ — U (I — 

2 n n 

G Atomic — -. — 1" -. — 77 (8) 

tU! + p tU! + p — U 
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FIG. 1: The zeroth moment of the self energy versus the density normalized to the exact value C/^ x nH^Hl^ This data is 
in 2-dimensions with U=10,W=2. The inset shows the k independence of the sum rule along the (11) direction for the case 
n=.04, with variations in the sixth significant figure. 
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FIG. 2: High frequency (UHB) DOS in 2D, U=10, n=l/20; As the hopping is decreased, the UHB feature does not become 
narrow, but rather maintains a width of 0(U). The sharp k-dependent features narrow as the hopping decreases. The broad 
continuum is essentially k-independent, showing very little dependence on the hopping in the limit of strong coupling. In this 
limit, the UHB becomes completely independent of the bandstructure. In figure[3]the broad UHB of the full band can be seen 
with the Hubbard- 1-like Gi superimposed. In Gi, the UHB feature is broadened only by rj and the LHB is suppressed for 
clarity. 



The breakup of the Greens function into two parts, with energies or ^ U and weights 1 — n/1 and n/2 is of course 
the fundamental factor that leads one to the picture of upper and lower Hubbard bands. Hubbard's contributior^i 
was to provide a Greens function for finite hopping W using an equation of motion method that extended the Atomic 
limit, although the details of his treatment came in for severe criticism- due to the failure of his scheme to ever yield 
a Fermi liquid with the Luttinger Ward^^ ordained Fermi surface. The present scheme of ladder diagrams achieves 
this interpolation smoothly and exactly, if only in the limit of low densities. From Fig. ([2]), we see that the sharp 
feature is accompanied by a broad background of width 0{U) that presumably arises from the uncontrolled O(n^) 
corrections to the ladder diagram self energy sum rule Eq. (jS]). 
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FIG. 3: 2D, U=10,W=2.5, n=l/20; The spectral function at three values of the wave vector (0,0), (^, ^), {it, it), in blue, red 
and gold colours. Besides the quasiparticles we observe three features emerging in each spectral function. Most obvious is the 
UHB feature which lies at a oj « 0{U) and integrates to a weight of n/2 + 0{n^). This feature is dramatically broadened in the 
self consistent G also becoming less k-dependent. On each edge of the quasiparticle band we observe small dispersing features. 
Ref. lfisl l have previously identified the negative frequency feature as a 2- hole antibound state while Ref.(f23) has discussed a 
particle-hole antibound state just above the quasiparticle band. These features are essentially unchanged in going from Gi to 
the exact G. 



C. Emergence and structure of the Split bands of Hubbard 



In the ladder diagrams, it is straightforward to identify the origin of the upper Hubbard band: the scattering 
amplitude r((3) at frequencies fig ~ U, has a pole in the first iteration, i.e. at the level of Gi with 

This pole was noted very early in works Ref. ([T2IIT3I ) who identified this pole as the origin of strong correlations and 
Gutzwiller type factors. In Fig. ([3]), we see that the spectral function obtained from the first iteration i.e. Gi shows 
a sharp feature at a higher energy of 0{U) that arises from this pole. This peak disperses and may be viewed as a 
"baby version" of the upper Hubbard band. Next a self consistent treatment of this theory with r{Q; [G]) evaluated 
with G (rather than Go) broadens the upper band substantially as seen in Fig. ([3]). It is interesting that the lowe 
Hubbard band, i.e. the structure at energies below U are stable with respect to the iterations, and are hardly different 
between the first iteration scheme and the final one. 

We also see in Fig. ([3]), the existence of two features that have been commented upon in literature. The feature near 
the band bottom that disperses, is the so called hole- hole bound state note by Randeria and Englebrecht Ref. (fisl ). 
whereas the hump near the leading edge is a particle hole bound state feature noted by Anderson Ref. (|23 ). These 
features coexist with the other, dominant ones, namely the quasiparticle peak of the Fermi liquid and the broadened 
upper Hubbard band peak. If we replace the log linear scale in Fig. with a linear linear scale as in Fig. the 
UHB becomes almost negligible compared to the LHB feature. 
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FIG. 4: 2D, U=.25, U=10, W=2, n=.049. The local pT,(i^) divided by U versus cj. The two chosen values of U are in the 
weak coupling (blue U = .25 ) and strong coupling (red U = 10) ranges respectively. We see at the lowest temperatures that 
the the self energy curves overlap when scaled by U displaying a characteristic quadratic dip at the chemical potential. 



D. Frequency Dependent Self energy 

We next display the self energy in the ladder scheme. The spectral density for the self energy is given in Eq. ([5]) , 
and it is possible to obtain an equation for its momentum sum, i.e. a local self energy density 

-^^Ps(fc,t^)= j dvpG,ioc{v)pr,ioc{i' + io){f{ixj) + nB{io + v)). (10) 

k 

For comparison, we note that the local self energy in the atomic limit considered in Section IIIBI is given by a single 
delta function centered at?7(l — ^) — /ias: 

PAtormciio) ^ -{1 - -) 5[u: + p - U [1 - -)]. (11) 

We also note the form of this object for a Fermi hquid at finite T 

PLocal (w) a W + SBackground(UJ), (12) 

a simple second order self consistent theory (corresponding to truncating the ladders at the first rung) gives the 
picture of this in a Fermi liquid Fig. (jj]). 

We see in Fig. [5] that the ladder scheme inherits both a quadratic minimum at w = from the Fermi liquid and a 
large and broad feature near uj ^ U from the emergent Hubbard upper band. The inset emphasizes the Fermi liquid 
aspect, and the reader will observe that the absolute scale of this function is dominated by the UHB feature. In Fig. 
|6]the density of states of the Greens function pQ{k,v) is illustrated, along with the real and imaginary parts of the 
self energy. The small feature in the DOS at the energy scale U is the UHB. We see that the real and imaginary parts 
of the self energy reflect its presence in a profound fashion, that would be hard to guess from the size of the peak. In 
detail, it is interesting that the real part of the self energy does display a linear behaviour in w with a known slope 
as one expects in the intermediate frequency range Q lo <^ U from the theory of extremely correlated electronic 
systems in Ref. (2426). 

When W = Q the UHB has a weight which is independent of momentum. However, for finite W , momenta near the 
top of the band will transfer weight more readily to the UHB. Fig. ([7]) illustrates this progression. We show in Fig. (|S]) 
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FIG. 5: The local self energy spectrum in 2D, U=10, W=2, n=.05 The log scale plot shows the full scale of the UHB. The 
inset highlights the quadratic minimum at low energies. The quadratic minimum drops below the scale of r] so it can be said 
to represent an infinite lifetime. 
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FIG. 6: The 2D DOS i.e. the momentum averaged spectral function paik, v) and the momentum averaged ps{k, v). The LHB 
feature is the sharp peak near a; ~ 0. The UHB feature in the DOS is nearly invisible here but lies just below the feature in 
Pe scaled down by a factor of oj^. The real part of the self energy for to >Q initially drops linearly with frequency over a range 
W ^ u) ^ as required in the limit of extreme correlations^!^. It then flips at the threshold of the UHB, rising across the 
range of the UHB until at the highest energy it begins to decay down towards the Hartree term at infinite energy. 
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FIG. 8: From the convolution structure of p^{k,uj) we see that the local objects of S and F are related by the ratio n/2 when 
LJ > W for all values of W/U. In the strong coupling limit where the upper band is essentially independent of k, this relationship 
will be approximately true for each wavevector. On the negative frequency side, the thermal function act differently such that 
the ratio for to < —W is approximately (1 + ^). 



that the behaviour of the local spectral function {pr{Q, v))q closely follows that of the local self energy pt.{v)- If we 
look at large uj such that we can make the approximation lo + v uj, the integral for p^{k^uj) in Eq. ([5]) reduces to 

■J^^P^ik,^) PT,loc{^^), (13) 



accounting for the similarity of these in Fig. |S1 
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FIG. 9: ID U=10 W=.56 (dashed), W=.196 (solid), n=.15 Id T=.005. Here mi is essentially the zero temperature quasiparticle 
occupation, while m2 accounts for the LHB particle addition spectrum. The sharp step in occupation occurs precisely at the 
Luttinger Fermi surface which satisfies the Luttinger Ward sum rule. The sum of mi and m2 is less than one due to the weight 
transferred to the upper band. The total lower band weight approaches l-n/2 as U/W goes to infinity. 



E. Momentum occupancy 

We next turn to the momentum occupancy nik = {c) {k)c{k))\ this can be obtained from the Greens function 
or pG{k,v) by integration over the frequencies. In order to understand and illustrate the nature of the LHB and 
UHB breakup of this important object, we carry out the integration up to the Hubbard-Mott gap energy ojg. This 
energy scale is well defined when W <^ U , and in case of smaller U W ii requires a definition. In our work, it is 
operationally defined as the energy where the spectral density {pdk, v))k is minimum. Thus we define three objects 
mj{k) with j = 1, 2, 3 

TOi(fc) / dujpG{k,uj) (14) 



m2{k) = / duipoikjU}) (15) 

/•oo 

m3(fc) = / dujpG{k,uj). (16) 



Here mi(k) represents the momentum space occupancy of the occupied states that lie below the chemical potential. 
These are automatically in the LHB for energetic reasons, and satisfy the sum rule rni(fc) — n/2 x Ng with a sum 
over the entire Brillouin zone (BZ). Next m2{k) represents the LHB contribution to the unoccupied states, since the 
chemical potential lies within the LHB. If we send [/ — ^ oo then we are left with only the LHB, and in that limit, 
we expect the sum mi(fc) + m2(fc) = 1 — § pointwise at each k. However for finite but large U this sum differs from 
1 — ^ by terms of 0{t/U), and the UHB comes into play. Indeed ra^{k) refers to precisely the UHB contribution to 
the momentum occupation, and its momentum average over the BZ is ^. These are displayed for typical parameters 
in Fig. The sum of all three m functions should add to unity for each wave vector. However, due to the finite 
frequency resolution of our numerics this sumrule is only approximately satisfied. We limit the error to < 1% by 
reducing our frequency step dio. The error is concentrated near kf where the spectral function is sharpest. 

In Fig. we display the k dependence of the three occupancy functions for a typical set of parameters. It is clear 
that the Luttinger Ward Fermi surface controls the variations of the functions mi and TO2, which complement each 
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FIG. 10: The doublon dynamics breaks into two regimes: a sharp decay at early times followed by a long exponential tail. The 
magnitude of the initial decay depends strongly on the density. In the limit n — the initial decay disappears, indicating that 
the UHB is comprised of sharp features only in the limit of vanishing density. In the right panel, the U dependence of the long 
time decay is shown, it slows down and is finally limited by the level broadening r; assumed in our numerics. 

other so that the sum is almost a constant. 



III. DOUBLONS AND THEIR DYNAMICS 
A. Doublon Decay in the low density limit 

In the recent experiments^'^ the Hfetime of doublons created by optical excitation of the trapped atoms has been 
carried out, providing us with an added impetus for this study. The experiments actually study the decay of a 
highly non equilibrium initial state \^initiai) with a finite fraction of excited doublons, i.e. {ipinitiai\D\^p initial) oc Ns, 
where the doublon number D = ni-\-nii. The object studied is the time evolution of such a state followed by a 
measurement of D and then a projection on to the evolved state i.e. 

= {ll^Initiall CXp {it^i/} D CXp {-it^iJ} I nitial) ■ (17) 

Here and below we use the symbol tr to denote real (Schrodinger) time, thus distinguishing it from the band hopping 
parameter t. Such a correlation function is not usually amenable to study near equilibrium type situations studied in 
many body physics. The initial state is itself quite far from being an equilibrium (ground) state. However, in the limit 
of very low densities, one can approximately view the initial state as the vaccuum or few particle state with a few 
doublon excitations- and within this picture we may ask how a single doublon decays. This is roughly the question 
of the lifetime of a state in the upper Hubbard band, and thus related to our general theme in this work. 

We are able to calculate the lifetime of a doublon within the ladder scheme, and hence presumably an exact answer 
at low densities as argued here. We next provide a discussion of the function 7 in a low density Fermi liquid. We 
start with the correlation function defined for Matsubara time r > in terms of the two particle Greens function^*' 

l{r, t) ^ G{l^^^(rr, rr; 0, 0) = (c,,t(r)c.,^(T)cS^^(0)cJ^^(0)), (18) 

and an analogous expression for real times 7(r, tr). This object can be expressed in terms of the scattering amplitude^l 

as 

7(r,i,) = ^ f dnpr{Q,i^){l + nB{i^))e-'Q^-"'*^. (19) 
Q •' 

In Fig. [TUl we display 7(0, t^) within the ladder scheme. As the density is increased, the UHB becomes broader 
and less k-dependent, however sharp k-dependent features persist with weight which decreases as n goes to zero. The 
k-dependent pieces remain sharp and determine the rate of the long time exponential decay. On the other hand 
the k-independent pieces, being broad, determine the short time decay. Due to our finite frequency resolution these 
numerics do not see the long time exponential decay becoming infinitely long once t < rj. 

We have also computed the off site correlation function 7(1, tr > 0), Fig. [TT] shows that even the site directly 
adjacent the created doublon has a very small amplitude. 
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B. Exact Solution of the Doublon Decay Problem for Two Particles. 

In addition to the discussion of the low density case, we are able to solve exactly the admittedly simple problem 
of the dynamics a single doublon in the Hubbard model, and from this study provide some feeling for the validity of 
the ladder scheme. The single doublon problem is solvable since for two particles of opposite spin, we have a total 
momentum quantum number and in each sector of this, we have a single particle type Schrodinger equation to solve. 
Let us first outline this problem and its solution with regard to the correlation function 

7(r, tr) = (0 I CrMU)CrM'^lliO)4iA^) I O)' (20) 

Here the average is with respect to the vaccuum state with no particles, although below we will use the average over 
the thermal distribution function for a low density Fermi liquid. In the case of two particles, it is in fact possible to 
show that 7(r, ir) is related to the correlator £,{tr) in Eq. (fT7|) exactly through 

atr) = Y.\^ir,U)\'- (21) 

r 

This follows upon using the fact that with only two particles in the system, the destruction operator Cr,t(ir)cr,4.(ir) 
can only connect to the vaccuum state. We expect this relation to be only approximately true for a dense Fermi 
system but useful since it can be computed with relative ease by one of several techniques. It is also dominated by 
the term r = as shown explicitly below in FigfTTl and hence it is useful to regard |7(0, tr)P as an estimator of £,{tr)- 
In Ref. (@), Demler et. al. estimate ^{0,tr) by an argument that is appropriate in an incoherent Fermi system, and 
estimate that this function decays on a time scale that is given as 

-=Atexp{-B (22) 

T W 

The vanishing of the rate as TV — )■ is expected in view of the conservation of the doublon number in the absence of 
electron hopping, the coefficients are estimated from experiments on the 3-d cubic lattice {W = 12t) as A '-^ .9 ± 0.5, 
and B - 1.6 ±0.16. 

For the two particle problem, we have exact analytical and numerical solutions. In the interesting case of U > W 
in d dimensional hypercubes with nearest neighbor hopping, we can write 

7(0, t,) = JL{0,U)+Ju{0,tr), 

^u{0,tr) - e-'(^+4'^V)*.-j„^(^t,), (23) 

where the LHB contribution 7^ ^ 0{{W/U)'^) and negligible. The second term arises from the UHB, and for 
intermediate W <^ U is related to the Bessel function Jg whereby it decays as a power law rather than as an 
exponential. This is understandable since the two body problem is an integrable system, and we expect that in the 
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FIG. 12: Doublon decay on a cubic lattice with [/ = 15 and W — 12. The shape of | ^u{Utr) P (red curve) initiaUy deviates 
shghtly from the exact numerical result (blue curve) due to the neglect of the 7l term, which decays much more quickly than 
the UHB contribution. 
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FIG. 13: Doublon decay on a cubic lattice with U = ^ and W — 12. In the case U < W, it is much more difficult to find an 
exact analytical form, so only the numerical result is displayed. 

low density limit, this power law would be replaced by an exponential type decay. The function can be found 
easily (see Appendix ) by numerical means and Figs ll2l and 1131 give us a picture of the decay. 

In Fig. [H we show that the Half Width at Half Max (HWHM) of the computed 7(0, ir-) leads to a rate ^^J^^^^ 
which has a behaviour that is similar to that in the experiments Eq. (j22p . 

IV. CONCLUSIONS 

In conclusion, we have shown that the self consistently computed ladder diagrams provide a detailed picture of the 
split bands for the Hubbard model. The UHB has a distinct shape that is captured here and related to the shape 
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1 2 3 4 5 

2zt/U 

FIG. 14: Two theoretical calculations, from ladder diagrams of Eq. (|19|) in 2-dimensions (blue) and the exact 2 particle solution 
from Eq. (|21[) and Eq. (|A10|) in 2- and 3- dimensions (red and gold). These are compared to the experiment Eq. (|22p in 
3-dimensions (green), scaled to coincide at weak coupling by a factor 26.4. The theory and expermient are in very different 
limits of physical parameters, but have a similar shape except at large U/t. 

of the two particle scattering amplitude. We have delineated how the lower Hubbard band occupation is influenced 
by the passage to large U. Here the background momentum occupance found in variational studies of the Gutzwiller 
approximationSS, arise here dynamically. Finally, we have shown that the decay of the doublon in such a system can 
be calculated by the ladder diagrams as well as by exact methods for very low densities, and the shapes of these 
curves are fairly close to those found in recent experiments on atomic traps performed under very different physical 
conditions. 



Acknowledgments 

This work is supported by DOE through a grant BES-DE-FG02-06ER46319. We are grateful to D. Huse, H. R. 
Kishnamurthy and M. Rigol for helpful discussions. 

Appendix A: Exact correlation functions for the two particle Hubbard Model 

We consider the Hubbard model with two particles, one spin up and the other spin down. Our goal is to calculate 
the following correlation function. 

l{tr) = (0 I C,^Qte-'^*'-cj^4 I 0) = lu{tr) + idtr)- (Al) 

The two parts arise from intermediate states that are in the two split bands. Thus 

lu{tr) = ^ |(H c\^cl I 0) p e-^-*- 

veUHB 

lL{tr) = E KHc.Vl, |0)pe-^-*'-. (A2) 

veLHB 

We now calculate the eigenvalues and eigenstates for the 2 particle Hubbard model. As our basis we take momentum 
eigenstates. 

|g,fc)^ct .t |0) (A3) 
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Q is the total momentum of the state, and both Q and k can be any vector in the first Brilluon zone. The Hamiltonian 
acts on the basis in the following way. 

H\Q,k) = Ek\Q,k)^^Y.\Q^P)^ 

p 

where — (eg-k + ik)- The Hamiltonian conserves total momentum. Thus, we can diagonalize each total momentum 
sector independently. Each sector will have Ns eigenstates, where Ng is the size of the lattice. We now fix Q and work 
in a particular total momentum sector. The basis states now depend on a single index k. The Ej- 's will in general be 
degenerate, and we take an E with degeneracy n, i.e. deg{E) = n, corresponding to states | Q,ki)... \ Q,kn)- From 
these we can make an n — 1 dimensional degenerate eigenspace of the Hamiltonian with energy E which we shall call 

I V')deg- 

n 

\il^)deg =^a^\Q,h) ^aj = (A5) 

i—l i 

One can see that these are eigenstates with energy E since potential energy term goes to zero due to the condition 
J2i cti = and the kinetic energy term gives E times the state. Suppose there are p unique values of E in this total 
momentum sector. 

deg{Ei) + ... + deg{Ep) = Ns (A6) 

By forming states in the way described above, we can obtain Ns—p eigenstates | ip)deg that are independent of U. We 
obtain the remaining non trivial (i.e. U dependent) p eigenstates by plugging the following state into the Hamiltonian. 

IV'q)=E*q(^)IQ'^) I Vq) =Aq I Vq). (A7) 

k 

Here we consider states with a fixed total momentum Q since this object is conserved. This yields the following results 

We can see explicitly from Eq. (jASp that {ipq \ 'ip)deg — since basis states with equal E have equal coefficients, and 
therefore the condition cti = makes this state orthogonal to the degenerate manifold of states in Eq. (jASp . There 
are p — 1 solutions of Eq. ()A8|) which lie in between the p distinct £"s. The corresponding states are in the lower 
Hubbard band. The | ip)deg found earlier also lie in the lower Hubbard band since these states are independent of U. 
There is one solution of Eq. (jA8l) for which Ag > E^ax and is of order U ii U > W. The corresponding state lies 
in the upper Hubbard band. Thus for each fixed Q sector, there is one state in the upper Hubbard band. We now 
consider the doublon state. 



Q,k 



We can rewrite 



where the Q in the above sum stands for the p states described by Eq. (|A8|) in the total momentum sector Q. Since 
{ipd I ip)deg — we didn't have to take the degenerate states into account when calculating the correlation function. 
Furthermore, we see that 

I (^Q I f= (All) 

where cq is from Eq. (jASp . 

QeUHB Q 
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In the above sum, each Q now represents only one state, since there is only one UHB state in each total momentum 
sector. We first evaluate this in one dimension, and then generalize to multiple dimensions. The sum can be turned 
into an integral. 

^uitr)^l r -j^e-^^^^'^dQ (A13) 
Converting Eq. (|A8p into integrals, we find that 



Aq = {U' + cos^ (A14) 



^([/^4-16t^cos^|-)5 (A15) 
For U > W, we keep corrections of 0{jp) in Aq and drop all corrections in Cq , yielding 

lu{tr) ^ - r e-*^(i+«^ ^^'-dQ (A16) 



2 4/^ 

lu{tr)^e~^(^+^^^''-Jo{—U) (A17) 

In two dimensions, Eq. (jA8l) becomes an elliptic integral so there is no closed form answer for the upper band 
eigenvalues in terms of elementary functions. However for U > W, keeping corrections to the same order as we did in 
deriving Eq. (|A16I) . we can easily generalize to higher dimensions. 

Ag^f/a + S^StiCos^ (A18) 

luiU) ~ e-(^+4'^^)*'- Jo^(^i.) (A20) 

The other contribution to 7(tr) is 7L(^r)- However, from degenerate perturbation theory, we know that provided 
U > W \ {ly \ ip)d P is 0{jj^) smaller for veLHB than it is for the upper Hubbard band. Hence, ^L{tr) is a small 
correction to ^u(tr)- 

l{tr)~lu{tr) (A21) 
I l{tr) P« Jn-^tr) (A22) 

In conclusion, the doublon decay in the 2 particle Hubbard model in the regime U > W is dominated by with 
the much faster decaying 7^ giving a small correction. To a good approximation, the shape of the decay of | 7(ir) P 

is jQ'''{^j-tr). 
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n(Q) 
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P1-P2 

where we used Eq. ([l}. The spectral density of 7 is now easy to find 

-pn 



Pt = 



(1 + t/i?en)2 + (7r(7pn)2 ' 
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